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Abstract 


The  use  of  finite-difference  and  finite-element  computer  codes  to  solve  problems  involving 
fast,  transient  loading  is  commonplace.  A  large  number  of  commercial  codes  exist  and  are 
applied  to  problems  ranging  from  fairly  low  to  extremely  high  damage  levels  (e.g.,  design  of 
containment  structures  to  mitigate  effects  of  industrial  accidents;  protection  of  buildings  and 
people  from  blast  and  impact  loading;  foreign-object  impact  damage;  design  of  space  structures 
to  withstand  impacts  of  small  particles  moving  at  hypervelocity,  a  case  where  pressures 
generated  exceed  the  material  strength  by  an  order  of  magnitude).  But,  what  happens  if  code 
predictions  do  not  correspond  with  reality?  This  report  discusses  various  factors  related  to 
material  interfaces  in  Lagrangian  and  Eulerian  shock  wave  propagation  codes  (hydrocodes), 
which  can  lead  to  disagreement  between  computations  and  experience.  Companion  reports  focus 
on  problems  associated  with  meshing  and  constitutive  models  and  the  use  of  material  data  at 
strain  rates  inappropriate  to  the  problem.  This  report  is  limited  to  problems  involving  fast, 
transient  loading,  which  can  be  addressed  by  commercial  finite-difference  and  finite-element 
codes. 

This  report  has  been  accepted  for  publication  in  a  future  volume  of  the  International  Journal 
of  Impact  Engineering. 
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1.  Introduction 


In  companion  papers  [1,2],  sources  of  error  in  calculations  involving  fast,  transient  loading  due 
to  various  meshing  options  available  in  commercial  wave  propagation  codes,  the  use  of  constitutive 
models,  and  use  of  data  at  strain  rates  inappropriate  to  the  problem  are  discussed.  In  this  report,  the 
effects  of  contact-impact  algorithms  in  Lagrangian  codes  and  the  presence  of  material  interfaces  in 
Eulerian  codes  are  discussed.  Attention  is  limited  to  conventional  production  Lagrangian  and 
Eulerian  hydrocodes  that  are  readily  available  in  industry  and  government  research  centers. 
Arbitrary  Lagrangian  Eulerian  (ALE)  codes  and  meshless  methods  (such  as  the  free  Lagrangian  or 
the  smooth  particle  hydrodynamics  [SPH]  technique)  are  not  covered  here. 

The  ability  to  treat  problems  where  adjacent  components  may  independently  slide,  separate,  and 
impact  along  free  surfaces  and  material  interfaces  is  crucial  in  many  technical  fields.  The  nuclear 
industry  analyzes  the  impact  of  shipping  casks  containing  radioactive  materials,  pipe-to-pipe 
impacts,  and  soil-structure  interaction  problems.  Crashworthiness  of  vehicles  is  another  important 
area,  as  is  foreign-object-impact  damage  in  aircraft.  Erosion  and  failure  of  materials  of  high-speed 
aircraft  flying  through  rain  and  hail  is  still  a  significant  problem.  Contact  impact  and  sliding  with 
and  without  friction  is  as  important  in  metal  forming  as  it  is  in  biomechanics.  In  weapons  design, 
relevant  problems  include  the  structural  response  of  gun-fired  projectiles,  various  bomb 
configurations,  and  a  variety  of  shaped  charge  designs. 

Spatial  discretization  using  either  finite-difference  or  finite-element  techniques  can  be  carried 
out  in  a  Lagrangian  or  Eulerian  framework  (Figure  1).  In  a  Lagrangian  framework,  the  grid  is 
embedded  in  the  material  and  distorts  with  it  so  that  the  computation  tracks  the  motion  of  elements 
of  mass.  In  the  Eulerian  approach,  the  grid  is  fixed  in  space  and  material  flows  through  it.  Eulerian 
codes  work  well  for  situations  where  large  distortions  are  encountered,  such  as  in  hypervelocity 
impact.  There  are  advantages  and  disadvantages  to  both  approaches.  Lagrangian  codes  are 
conceptually  straightforward.  Since  the  grid  distorts  with  the  material,  time  histories  are  easily 
obtained  and  material  interfaces  are  sharply  defined.  However,  considerable  mesh  distortion, 
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Figure  1.  (a)  Lagrangian  Computational  Grid  and  (b)  Euierian  Computational  Grid. 
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entanglement,  and  nodal  intrusions  across  slide  lines  can  occur  and  special  logic  is  required  for 
modeling  deep  penetration  problems.  In  Eulerian  codes,  material  is  transported  from  cell  to  cell  in 
the  fixed  grid.  Considerable  diffusion  can  occur,  and  special  logic  is  required  if  interfaces  are  to  be 
defined  with  greater  precision  than  one  or  two  cell  dimensions.  Additional  information  can  be  found 
in  Zukas  et  al.  [3],  Zukas  [4],  Walters  and  Zukas  [5],  Hallquist  et  al.  [6],  Belytschko  and  Hughes  [7], 
Schwer  et  al.  [8],  and  Kulak  and  Schwer  [9]. 

2.  Lagrangian  Calculations:  Influence  of  Contact  Surfaces 
on  Solutions 

Figure  2  [  1 0]  shows  a  typical  interface  setup  between  two  contacting  bodies,  one  side  designated 
as  master,  the  other  as  slave.  Sliding-interface  algorithms  can  be  divided  into  three  broad  categories: 

(1)  the  penalty  method,  or  restoring  force  approach,  particularly  popular  in  structural 
dynamics  and  low-velocity  applications,  where  nodal  intrusions  are  corrected  over 
several  computational  cycles; 

(2)  the  “put-back”  logic  approach  used  in  most  hydrocodes,  with  or  without  erosion 
algorithms,  where  intruding  nodes  are  repositioned  on  the  contact  surface  in  a  single 
computational  cycle;  and 

(3)  novel  approaches,  such  as  the  pinball  algorithm  [11]. 

In  addition  to  these,  there  have  been  various  ad  hoc  approaches  used  for  specific  problems.  Each 
of  these  has  proven  useful  in  a  specific  spectrum  of  nodal  (therefore  impact)  velocities,  but  none  is 
sufficiently  general  to  work  in  all  applications  from  elastic  impact  through  hypervelocity  impact, 
where  perforations  occur  and  large  debris  clouds  are  formed. 

At  the  outset  of  a  Lagrangian  calculation,  where  interpenetration  of  two  or  more  bodies  can  be 
expected,  the  surface  nodes  likely  to  make  contact  are  designated  as  slave  nodes  on  one  of  the  bodies 
and  master  nodes  on  another.  With  modem  contact-impact  algorithms,  the  designation  is  more  or 
less  arbitrary.  These  master  and  slave  nodes  must  be  designated  by  the  code  user.  Codes  such  as 
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Figure  2.  Typical  Lagrangian  Contact  Interface  Showing  Master-  and  Slave-Node 
Designations  [10]. 
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ZeuS  [  1 2]  eliminate  this  nuisance  by  automatically  designating  all  geometric  boundaries  and  material 
interfaces  as  contact  surfaces.  The  user  can  remove  them  if  they  are  not  needed  but  is  never  required 
to  specify  them.  The  goal  of  sliding  interface  logic  (also  known  as  “contact-impact  algorithms”)  is 
to  prevent  intrusion  of  a  slave  node  through  a  master  surface.  Intrusions  must  be  corrected.  The 
procedures  are  outlined  in  Tables  1  and  2. 

The  most  comprehensive  contact  algorithms  of  this  sort  are  in  the  DYSMAS-L  code  [13]. 
Following  the  scientific  method  to  its  logical  end,  the  authors  of  the  code  exhaustively  considered 
every  possibility  of  slide-line  intrusion  and  programmed  against  it.  Because  all  potential  types  of 
intrusion  are  checked,  considerable  time  is  expended  in  DYSMAS-L’s  contact  processor  at  each 
cycle,  making  for  a  slow-running  code,  but  one  without  any  significant  errors  due  to  the  contact 
processor.  In  the  United  States,  the  approach  in  codes  such  as  DYNA  [14,  15],  EPIC  [16,  17], 
PRONTO  [18],  ZeuS,  and  others  has  been  more  pragmatic.  These  guard  against  the  major  sources 
of  error  that  are  likely  to  occur.  Once  intrusion  of  the  master  surface  is  detected,  the  offending  node 
is  repositioned  in  one  or  another  fashion  and  contact  conditions  and  momentum  balance  are  enforced 
locally  on  the  sliding  interface.  Pathological  cases  are  then  treated  as  they  occur. 

By  constrast,  with  penalty  methods,  the  restoring  force  on  an  intruded  slave  node  is  very 
sensitive  to  the  stiffness  coefficient,  which,  for  quite  some  time  in  the  history  of  Lagrangian  code 
development,  was  a  user-defined  parameter  and  a  rather  significant  source  of  error  if  the  code  user 
was  not  familiar  with  the  physical  problem  and  his/her  code.  The  current  practice  is  to  determine 
the  restoring  force  modulus  from  the  properties  and  geometry  of  the  element  hosting  the  sliding 
segment.  Prior  to  this,  results  could  be  dramatically  affected  by  a  poor  choice  of  stiffness  constant 
and  that  choice  was  very  much  a  matter  of  experience.  By  way  of  example,  Goudreau  and 
Hallquist  [19]  incorporated  procedures  in  NIKE  and  DYNA3D  to  compute  a  unique  modulus  for 
each  slave  and  master  segment  based  on  the  thickness  and  bulk  modulus  of  the  element  in  which  it 
resides.  They  state  that  . .  in  our  opinion,  the  lack  of  user  control  over  this  very  critical  parameter 
greatly  increases  the  success  of  the  method.” 
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Table  1.  Sliding  Interface  Approach  in  Lagrangian  Systems 


Designate  master  and  slave  material. _ 

Decompose  u  ,  u  into  normal  and  tangential  components. 

In  normal  direction: 

-motion  coupled  during  contact,  and 
-independent  when  separated. 

In  tangential  direction: 

-independent  when  materials  separated  or  interface  frictionless,  and 
-modified  when  there  is  contact  or  frictional  force  acts. _ 

Tied  sliding  -  maintain  displacement  compatibility  at  abrupt  grid  changes. 


Table  2.  Computational  Procedure  for  Sliding  Interfaces 


(1)  Identify  a  series  of  nodes  that  make  up  the  master  surface. 

(2)  Identify  a  series  of  nodes  that  make  up  the  slave  surface. 


uations  of  motion  to  both  master  and  slave  nodes. 


(4)  Check  for  interference  between  slave  nodes  and  master  surface: 
-define  search  region, 

-check  each  slave  node  not  on  master  surface,  and 

-if  penetration  occurs,  move  slave  node  to  master  surface. 


(5)  Once  slave  nodes  are  moved  onto  master  surface: 

-invoke  momentum  balance, 

-apply  frictional  forces,  and 

-open  voids  (if  tensile  forces  present). _ 

(6)  Repeat  (3>— (5)  for  each  At 


While  this  restorative  procedure  of  the  penalty  method  gave  a  smoothness  to  the  calculation  due 
to  the  fact  that  restoration  of  intruding  nodes,  and  therefore  element  volume  changes,  was  gradual, 
the  procedure  did  not  lend  itself  to  impacts  at  ordnance  velocity  or  hypervelocity  because  intrusions 
could  exceed  element  dimensions.  Hence,  most  existing  Lagrangian  hydrocodes  use  the  put-back 
logic  mentioned  previously.  The  intruding  node  or  nodes  are  put  back  onto  the  master  surface  in  one 
time  step,  and  momentum  balance  is  invoked.  Efficient  search  procedures  for  slave-node  intrusion 
are  required  since  a  Lagrangian  code  is  made  or  broken  by  the  accuracy  and  efficiency  of  its  contact 
logic. 

Ideally,  the  operation  of  the  contact  processor  should  be  transparent  to  the  user.  Very  much, 
however,  depends  on  the  gridding  in  the  vicinity  of  the  contact  surface  and  the  local  time  step. 
Breidenbach  et  al.  [20]  provide  a  comprehensive  discussion  of  pathological  cases  that  can  arise  with 
slide-line  logic.  Grove  et  al.  [21],  studying  tungsten  projectiles  penetrating  titanium  targets, 
exercised  the  various  erosion  options  in  EPIC  (upon  failure  of  an  element,  its  connectivity  is  no 
longer  considered  and  the  nodal  masses  associated  with  that  element  can  be  retained  in  the 
calculation  or  ignored)  and  compared  results  using  the  data  for  the  Johnson-Cook  models  [22, 23] 
in  the  EPIC  materials  library  with  results  using  strength  data  reported  by  Burkins  et  al.  [24] .  Results 
for  one  of  the  calculation  sets  are  shown  in  Figure  3.  The  differences  between  the  calculation  sets 
are  due  to  a  combination  of  data  for  the  selected  constitutive  model  and  the  behavior  of  the  sliding 
interface  logic  when  dealing  with  free-flying  masses. 

The  Lagrangian  codes,  which  first  appeared  in  the  early  1 960s,  such  as  HEMP  and  TOODY,  had 
no  means  of  dealing  with  the  large  distortions  that  can  occur  under  impact  and  explosive  loading. 
They  required  that  the  contact  surface,  or  sliding  interface,  specified  at  the  beginning  of  a  calculation 
remain  unchanged.  In  effect,  an  impenetrable  membrane  was  placed  over  and  attached  to  the 
colliding  bodies.  Thus,  colliding  bodies  could  not  fail  but  could  distort  ad  infinitum.  A  dramatic 
example  of  such  distortions  is  shown  in  Figure  4  [25].  Thus,  these  codes  were  limited  to  the  study 
of  the  early  stages  of  impact  phenomena.  In  order  to  extend  the  capability  of  these  codes,  rezoning 
techniques  were  invented.  In  rezoning,  a  new  Lagrangian  computational  mesh  is  overlaid  on  the 
distorted  mesh.  A  rezone  algorithm  maps  mesh  quantities  of  the  severely  distorted  mesh  onto  the 
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Figure  3.  Comparison  of  Eroding  Slide-Line  Options  on  Numerical  Results  for  a 
Deep-Penetration  Problem.  The  Type  =  - 1  Option  Retains  the  Mass  of  the  Eroded 
Elements  at  the  Nodes,  While  the  Type  =  -2  Option  Discards  Mass  Once  an 
Element  Erodes.  Also  Compared  Are  Simulations  With  Strength  Data  Taken  From 
the  EPIC  Library  to  Simulations  Using  Strength  Data  Reported  Elsewhere.  The 
Solid  Line  Represents  the  Experimentally  Obtained  Data  [21]. 
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URGE  MESH  DISTORTIONS— URGE  TRUNCATION  ERRORS. 


Figure  4.  Example  of  Extreme  Distortion  in  a  Lagrangian  Calculation  [25]. 
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new  mesh  such  that  conservation  of  mass,  momentum,  total  energy,  and  the  constitutive  relationship 
are  satisfied.  Rezoning  was,  and  remains,  a  nontrivial  process.  For  deep  penetration  problems,  30 
to  50  rezones  are  not  uncommon.  Since,  at  each  rezone,  the  accumulated  history  of  several  elements 
is  combined  into  a  larger  element,  frequent  rezoning  renders  the  computational  mesh  semi-Eulerian 
in  that  large  distortions  are  realized  but  material  history  and  location  of  material  boundaries  are 
diffused.  Indeed,  if  rezoning  were  to  take  place  at  each  cycle,  the  computations  could  be  effectively 
Eulerian.  Frequent  rezoning  gives  computational  results  a  smoothness  characteristic  of  Eulerian 
calculations.  Many  codes  today,  such  as  DYNA,  incorporate  automatic  rezoners.  Despite  the  label, 
they  still  require  occasional  human  intervention.  This  suggests  that,  despite  good  intentions,  rezoned 
calculations  can  be  made  to  mirror  the  expectations  of  the  user  monitoring  the  process.  An  example 
of  a  rezoned  calculation  is  shown  in  Figure  5  [26]. 

The  eroding  slide-line  concept  proved  to  be  a  major  breakthrough  in  extending  the  capabilities 
of  Lagrangian  codes  to  do  problems  involving  deep  penetration  or  multiple-plate  targets  and  other 
situations  where  damage  is  highly  localized.  The  term  “erosion”  in  this  context  does  not  refer  to  a 
physical  failure  mechanism.  Rather,  it  describes  bookkeeping  procedures  that  allow  dynamic 
redefinition  of  master  and  slave  surfaces  when  very  large  distortions  occur.  The  overall  result  makes 
it  appear  that  the  colliding  bodies  are  “eroding.”  Some  examples  of  the  success  achieved  with  this 
technique  are  shown  in  Figures  6-8.  Figure  6,  taken  from  Kimsey  and  Zukas  [27],  shows  the 
computational  grid  overlaid  on  the  recovered  target  from  the  impact  of  a  long  steel  rod  onto  a  semi¬ 
infinite  armor  steel  target  at  a  velocity  of  3.1 14  km/s.  Figure  7  [28]  shows  debris-cloud  formation, 
propagation,  and  interaction  with  other  components  of  a  space  vehicle.  Figure  8  [29]  shows  the 
penetration  of  a  50-ply  aluminum  plate  by  a  steel  projectile  15  ps  after  impact. 

Eroding  slide-line  techniques  are  available  in  virtually  all  Lagrangian  hydrocodes  today.  The 
name  of  the  procedure  may  change  (e.g.,  “advected  elements”  is  popular  today),  but  the  procedure 
is  essentially  as  described  previously.  While  very  useful,  such  techniques  should  be  used  with 
caution.  Energy  is  not  conserved  in  most  sliding-interface  prescriptions.  The  total  energy  in  each 
calculation  should  therefore  be  carefully  monitored.  Losses  of  4-10%  in  contact  calculations  can 
be  tolerated  for  most  high-velocity  impact  problems.  Anything  greater  than  10%  is  usually 
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(a)  Mesh  Before  Rezoning 


(b)  Mesh  After  Rezoning 


Figure  5.  Lagrangian  Calculation  Before  and  After  Rezoning  [26]. 


Figure  6.  Overlay  of  Computed  and  Experimental  Hole  Profiles  for  Steel-Steel  Penetration, 
Vs  -  3.114  km/s  [27]. 
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Figure  7.  Damage  Due  to  Sphere  Impact  at  30°  [28]. 
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Figure  8.  Penetration  of  a  50-Ply  Aluminum  Target  15  ps  After  Impact  [29]. 
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indicative  of  a  serious  problem  and  should  prompt  a  check  of  both  the  computational  model  and  the 
data  used  for  the  constitutive  models. 

In  penetration  problems,  a  situation  often  encountered  is  the  continuation  of  penetration  after 
the  projectile  has  eroded.  This  is  shown  in  Figures  9  and  10.  So  many  nodes  accumulate  at  the 
centerline  element  of  the  target  (close  to  180  nodes  in  this  case)  that  the  element  can  no  longer 
withstand  the  impacts  and  fails.  As  it  fails,  more  nodes  fall  into  the  crater  and  the  target  appears  to 
unzip  along  its  centerline.  Above  this  artificial  crack,  the  crater  dimensions  and  penetration  depth 
are  about  as  expected.  The  cause  for  this  behavior  is  unknown  at  present,  but  clues  to  the  problem 
may  be  found  in  the  work  of  Brach  [30]  and  Walsh  [31].  Hydrocodes  are  formulated  assuming 
continuum  behavior.  The  impact  of  single  or  multiple  mass  points  against  an  intact  element  may 
violate  that  assumption.  Brach  gives  computational  procedures  for  discrete  particle  impacts  against 
continuum  structures.  Walsh  studies  the  effects  of  the  impact  of  minuscule  mass  points  against 
continuum  targets  and  found  that  the  apparent  strength  of  the  target  increased  by  a  factor  of  4.7  over 
pretest-measured  strength  properties.  The  formulation  of  the  slide-line  algorithm  may  also  be  a 
factor.  Introducing  an  analytical  formulation  for  the  contact  processor  in  version  4.23  of  ZeuS  (the 
calculations  shown  here  are  with  the  standard  version  4.22)  helped  to  considerably  reduce  the 
unzipping  effect. 

3.  Eulerian  Calculations:  Effects  of  Material  Interface  on 
Solutions 

In  Lagrangian  codes,  the  mesh  moves  with  the  material,  whereas  in  Eulerian  codes  the  materials 
flow  though  the  mesh.  This  is  typically  done  in  two  phases.  In  the  first  phase,  the  mesh  is  allowed 
to  distort  as  the  problem  is  advanced  in  time  (the  Lagrangian  phase),  and,  in  the  second  phase,  the 
distorted  mesh  is  remapped  back  to  the  original  mesh  (the  advection  phase).  Operator  splitting  is 
typically  employed  during  the  advection  phase  such  that  the  remapping  occurs  one  coordinate 
direction  at  a  time,  and  the  order  is  reversed  each  cycle.  For  example,  the  order  in  which  advection 
occurs  may  be  in  the  x,  y,  then  z  directions  one  cycle  and  in  reverse  order  the  next.  Similar  to  the 
contact  algorithm  in  Lagrangian  codes,  the  advection  phase  is  the  most  computationally  intensive 
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Figure  10.  Residual  Penetration  by  Eroded  Mass  Points. 
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in  Eulerian  codes.  Some  codes  (such  as  MESA  [32])  allow  subcycling  of  the  Lagrangian  phase  (i.e., 
allowing  several  Lagrangian  phases  before  performing  a  single  advection  phase). 

As  is  shown  in  Figure  1(b),  material  interfaces  are  not  well  defined  in  Eulerian  codes.  The 
interface  between  two  dissimilar  materials  is  known  at  best  to  within  some  fraction  of  a  cell 
dimension.  These  interface  cells  contain  materials  of  both  the  bodies  involved  in  contact  and  are 
therefore  designated  as  “mixed  cells.”  Not  knowing  the  location  or  orientation  of  an  interface  within 
a  mixed  cell  is  a  serious  drawback,  and  a  number  of  techniques  have  been  invented  to  overcome  this 
problem. 

The  earliest  attempt  to  refine  the  location  of  a  material  interface  in  an  Eulerian  calculation 
involved  the  use  of  massless  Lagrangian  tracer  particles.  A  number  of  these  would  be  introduced 
at  the  outset  of  the  calculation,  typically  defining  the  boundaries  of  colliding  objects.  This  required 
that  both  Lagrangian  and  Eulerian  calculations  be  performed  each  cycle.  The  massless  tracers  were 
connected  by  straight  lines.  This  approach  worked  initially,  but  problems  arose  as  the  distance 
between  tracers  increased  with  increasing  deformation.  Additional  tracers  would  need  to  be 
introduced  by  the  code  user  at  various  restart  times  to  prevent  materials  from  flowing  across 
boundaries.  In  effect,  this  amounted  to  manually  assisted  rezoning  in  codes  such  as  HELP  [33].  As 
an  example  of  this  approach,  a  typical  HELP  code  calculation  is  shown  in  Figure  11.  While  the 
addition  of  tracers  alleviated  the  problem  of  boundary  transgression,  it  also  introduced  a  subtle  error 
in  the  calculation.  The  placement  of  the  tracers  influenced  the  further  computation  of  displacements 
so  that,  all  too  often,  the  final  resultant  deformation  was  a  reflection  of  the  preconception  of  the  code 
user.  If  the  user  was  inexperienced  and  unfamiliar  with  high-velocity  impact  physics,  the  errors 
introduced  into  the  calculation  could  be  quite  far  from  subtle.  In  that  era,  there  was  a  very  small 
body  of  experienced  computational  analysts  who  performed  such  calculations.  Novices  tended  to 
serve  lengthy  apprenticeships  under  experienced  code  users,  many  of  whom  had  studied  under  code 
developers.  Hence,  despite  slow  computers  and  long  computation  times,  results  tended  to  be  quite 
good. 
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Figure  11.  A  Typical  Plane-Strain  HELP  Code  Simulation  Compared  With  Experiment. 
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Considerable  improvements  in  the  use  of  Lagrangian  particles  to  define  interfaces  were  made 
in  the  SMITE  code  [34-37].  In  SMITE,  a  second-order  Eulerian  finite-difference  code,  each 
material  had  its  own  independent  grid.  Thus,  the  mesh  spacing  and  number  of  mesh  grids  in  one 
material  is  not  affected  by  any  other  material.  The  extent  of  each  domain  was  determined  by 
particles  that  defined  the  domain  boundary.  These  points  were  moved  in  a  Lagrangian  sense  by 
integrating  the  ordinary  differential  equations  relating  their  positions  or  velocities.  Additional 
particles  were  added  as  needed  to  ensure  uniform  spacing  throughout  the  calculation.  The 
interaction  between  materials  was  through  the  boundaries.  The  boundary  points  were  subject  to 
ffee-surface  and  interface  conditions  and  provided  the  only  communication  between  the  various 
material  domains. 

The  next  step  in  complexity,  if  not  necessarily  accuracy,  involves  the  use  of  interface 
reconstruction  algorithms.  One  of  the  most  successful  (because  of  its  simplicity)  material  interface 
reconstruction  algorithms  was  the  first-order  accurate  simple  line  interface  calculation  (SLIC) 
algorithm  of  N oh  and  W oodward  [38].  The  algorithm  treats  each  coordinate  direction  independently 
and  represents  material  interfaces  as  straight  lines  either  perpendicular  or  parallel  to  a  coordinate 
direction  based  on  the  volume  fractions  of  neighboring  cells.  The  order  in  which  each  coordinate 
direction  is  treated  is  reversed  each  computational  cycle.  The  SLIC  algorithm  did  well  when 
material  flow  was  primarily  in  a  coordinate  direction  but  tended  to  reorient  the  material  interfaces 
if  material  flow  was  not  parallel  to  a  coordinate  direction.  For  example,  SLIC  would  tend  to  reorient 
the  material  interfaces  of  a  circle  moving  the  mesh  at  45°  such  that  it  would  take  on  a  diamond  shape. 
An  example  of  the  algorithm  is  provided  in  the  “balls-and-jacks”  problem  shown  in  Figure  12  in 
which  two  balls  (resembling  a  “figure  eight”)  and  two  jacks  are  moved  through  the  mesh  at  45  °  with 
x  and  y  velocity  components  of  100  km/s. 

Most  modem  production  Eulerian  hydrocodes  use  a  second-order-accurate  material  interface 
reconstruction  scheme.  One  such  scheme  is  Youngs’  interface  reconstruction  algorithm  [39].  The 
Y oungs  ’  algorithm  also  uses  lines  (or  planes  in  three  dimensions)  to  represent  the  material  interfaces; 
however,  the  lines  can  be  at  any  angle.  A  problem  with  the  Youngs’  algorithm  is  that  the  code  user 
is  required  to  input  the  order  in  which  materials  will  advect.  The  wrong  choice  of  material  advection 
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Figure  12.  The  Balls-and- Jacks  Problem  Using  the  SLIC  Algorithm,  (a)  Material  Interface 
Plot  Shows  the  Initial  and  Final  Configurations  of  the  Balls  and  Jacks.  The 
Interfaces  Are  Reoriented  in  the  Direction  of  Motion,  (b)  Material  Volume 
Fraction  Plot  Shows  That  Small  Fragments  of  Material  Are  Left  Behind. 
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order  can  lead  to  strange  results,  as  shown  in  the  balls-and-jacks  problem  in  Figure  13.  In  very 
complex  problems  involving  many  materials,  considerable  skill  is  needed  in  choosing  the  material 
advection  order.  Additionally,  the  advection  order  of  materials  in  a  complex  problem  is  rarely 
uniform  over  the  entire  computational  domain.  The  problems  with  the  choice  of  material  advection 
order  were  largely  solved  with  the  Sandia  modified  Youngs’  reconstruction  algorithm  (SMYRA) 
[40]  used  in  the  CTH  code  [41],  which  is  essentially  Youngs’  algorithm,  but  material  volume 
fractions  of  neighboring  cells  are  used  to  determine  the  advection  order  in  mixed  cells.  An  example 
of  the  balls-and-jacks  problem  using  SMYRA  is  shown  in  Figure  1 4.  A  more  comprehensive  review 
of  material  reconstruction  algorithms  can  be  found  in  the  work  by  Benson  [42]. 

Even  with  a  perfect  material  reconstruction  algorithm,  problems  can  occur  when  objects  with 
the  same  material  identity  come  into  contact.  The  interface  in  this  case  will  simply  disappear,  and 
the  material  behaves  as  though  it  were  bonded  together.  This  occurs  because  material  interfaces  can 
only  be  determined  between  different  materials  or  materials  and  void,  not  like  materials  (meaning 
materials  that  are  assigned  the  same  material  number).  An  example  is  shown  in  Figure  1 5  for  a  deep 
penetration  problem  involving  a  tungsten  alloy  projectile  penetrating  a  steel  target.  In  this  case,  the 
erosion  products  of  the  long-rod  penetrator  came  in  contact  with  the  penetrator  and  the  interface 
between  the  penetrator  and  erosion  products  can  no  longer  be  distinguished.  The  end  result  was  that 
the  simulation  underpredicted  final  penetration  due  to  the  penetrator  being  decelerated  from  the 
contact. 

Velocities  in  most  Eulerian  hydrocodes  are  either  face  centered  or  node  centered,  while  other 
cell  quantities  are  cell  centered.  Therefore,  materials  in  a  mixed  cell  all  have  the  same  velocity  field. 
Thus,  a  cell  moves  as  a  distorting  block  implying  a  no-slip  condition.  Therefore,  Eulerian  codes  do 
not  handle  problems  very  well  where  sliding  between  materials  occurs  or  where  materials  separate. 
Unfortunately,  sliding  occurs  in  most  problems.  Attempts  have  been  made  to  correct  this 
shortcoming.  For  example,  Silling  [43]  developed  a  boundary  layer  algorithm  for  sliding  interfaces 
(BLINT)  for  two-dimensional  problems.  The  algorithm  creates  a  “slip  layer”  in  a  user-defined  “soft” 
material.  The  strength  of  the  cells  in  the  slip  layer  is  set  to  zero,  allowing  sliding  to  occur  outside 
of  the  mixed  cells.  The  model  has  been  used  with  some  success  in  modeling 
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Figure  13.  The  Balls-and-Jacks  Problem  Using  Y oungs’  Interface  Reconstruction  Algorithm 
Using  Improperly  Chosen  Material  Advection  Orders,  (a)  Material  Interface  Plot 
Shows  That  the  Algorithm  Does  a  Better  Job  Keeping  Track  of  the  Material 
Interfaces  Than  the  SLIC  Algorithm  of  Figure  12(a).  (b)  Material  Volume 
Fraction  Plot  Reveals  That  Material  Incorrectly  Advected  Is  Being  Left  Behind. 
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Figure  14.  The  Balls-and-Jacks  Problem  Using  SMYRA.  (a)  Material  Interface  Plot  Shows 
Results  Similar  to  the  Youngs’  Algorithm  of  Figure  13(a).  (b)  Material  Volume 
Fraction  Plot  Reveals  That  No  Material  Fragments  Are  Being  Left  Behind. 
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Figure  15.  Example  of  a  Deep  Penetration  Simulation  in  Which  the  Erosion  Products  of  the 
Rod  Come  Into  Contact  With  the  Penetrator.  At  the  Point  of  Contact,  an  Interface 
Is  Indistinguishable.  The  Contact  Caused  the  Penetrator  to  Decelerate  More 
Rapidly  Than  Usual,  Causing  the  Simulation  to  Underpredict  Penetration. 
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rigid-body  penetrations  [44—46]  and  has  recently  been  extended  to  three  dimensions.  An  example 
problem  comparing  simulations  with  and  without  the  BLINT  model  is  shown  in  Figure  16.  Walker 
and  Anderson  [47]  attempted  to  overcome  the  problem  by  developing  a  model  for  multimaterial 
velocities  in  mixed  cells.  The  model  defined  a  cell-centered  velocity  and  allowed  each  material  to 
have  its  own  velocity,  which  was  advected  as  a  state  variable.  The  model  was  strictly  geometrical. 
The  authors  attempted  to  model  a  rigid-body  penetrator  problem  using  the  model  with  only  limited 
success. 

Because  the  cell  stress  increments  are  calculated  from  velocity  gradients,  using  the  face-centered 
or  node-centered  velocities,  a  single  yield  strength  must  be  calculated  for  mixed  cells.  Some 
Eulerian  hydrocodes  allow  the  user  to  select  the  formulation  for  determining  the  yield  strength  of  a 
mixed  cell.  As  an  example,  the  CTH  hydrocode  gives  the  user  three  choices  on  how  the  yield 
strength  of  a  mixed  cell  is  determined.  They  include  (1)  setting  the  strength  of  a  mixed  cell  to  zero; 

(2)  multiplying  the  yield  strength  of  each  material  by  its  volume  fraction  and  adding  them;  and 

(3)  multiplying  the  yield  strength  of  each  material  by  its  volume  fraction,  summing  them  and 
dividing  by  the  total  volume  fraction  of  all  materials  within  the  cell.  All  these  methods  for 
determining  the  yield  strength  in  a  mixed  cell  can  lead  to  calculating  unrealistic  results  when  the 
materials  within  a  mixed  cell  have  vastly  different  strengths.  Figure  17  compares  simulations  using 
the  previously  mentioned  strength  treatments  for  mixed  cells.  For  each  strength  formulation,  the 
penetrator  did  not  perforate  the  target,  when,  experimentally,  it  was  determined  it  should.  The 
reason  is  partly  due  to  the  strength  treatment  artificially  weakening  the  penetrator  material  in  mixed 
cells  and  partly  due  to  the  fracture  model  chosen.  If  a  solid  were  in  a  mixed  cell  with  a  gas,  the  solid 
can  undergo  unrealistic  deformation.  The  strength  treatment  used  in  mixed  cells  is  one  reason  most 
Eulerian  codes  cannot  model  rigid-body  penetration. 

Mixed  cells  can  also  cause  thermodynamic  problems.  This  occurs  because  most  Eulerian  codes 
typically  determine  the  thermodynamic  state  of  the  mixed  cell  on  a  subcell  level.  The 
thermodynamic  state  of  neighboring  cells  is  not  used  to  determine  the  state  of  materials  within  a 
mixed  cell.  Assumptions  must  therefore  be  made  to  give  a  closed  set  of  equations  to  determine  the 
state  within  a  mixed  cell.  The  assumptions  most  often  used  include  (1)  all  materials  are  at  the  same 
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(b) 

Figure  16.  Comparison  of  a  Deep  Penetration  Simulation  With  and  Without  the  BLINT 
Model,  (a)  BLINT  Model  Simulation  Shows  the  Ogival-Nose  Rod  Penetrating  as 
a  Rigid  Body.  The  White  Layer  Around  the  Rod  Shows  the  Slip  Layer  Where 
Sliding  Occurs,  (b)  Simulation  Without  the  BLINT  Model  Incorrectly  Shows  the 
Ogival-Nose  Rod  Eroding. 
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Figure  17.  Comparison  of  Simulations  Using  the  Different  Strength  Options  in  Mixed  Cells. 

All  Simulations  Show  Excess  Deformation  of  the  Rod  and  Incorrectly  Predict  That 
the  Penetrator  Will  Not  Perforate  the  Target,  (a)  Yield  Strength  of  Mixed  Cell 
Was  Determined  by  Multiplying  the  Volume  Fractions  of  Materials  That  Can 
Support  Shear  With  the  Materials’  Yield  Strength  and  Summing  Them,  (b)  Yield 
Strength  of  Mixed  Cells  Was  Determined  by  Multiplying  the  Volume  Fractions  of 
Materials  That  Can  Support  Shear  With  the  Materials’  Yield  Strength,  Summing 
Them,  and  Dividing  by  the  Sum  of  Volume  Fractions  of  Materials  Supporting 
Shear,  (c)  Yield  Strength  in  Mixed  Cells  Is  Zero. 
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pressure  and  temperature;  (2)  all  materials  are  at  the  same  pressure,  but  the  temperature  of  each 
material  may  vary;  or  (3)  all  materials  may  have  their  own  pressure  and  temperature.  Assuming  that 
all  materials  have  the  same  pressure  as  in  the  first  two  assumptions  is  obviously  incorrect.  While 
pressures  may  eventually  equilibrate,  this  usually  happens  over  time,  not  instantaneously. 
Additionally,  assuming  all  materials  have  the  same  pressure  can  cause  problems  to  arise  when 
fracture  is  to  be  modeled  using  a  tension  criteria  in  a  cell  containing  a  solid  and  fluid.  The  solid  will 
never  fail  because  fluids  do  not  support  tension.  Assuming  thermal  equilibrium  implies  an  infinite 
thermal  conductivity.  This  assumption  was  used  in  the  CSQII  code  [48].  The  assumption  that  all 
materials  have  their  own  temperature  and  pressure  on  the  surface  is  the  most  correct;  however, 
additional  approximations  have  to  be  made.  These  include  assuming  the  volumetric  strain  is 
constant  among  all  materials  in  a  mixed  cell.  This  means,  in  a  mixed  cell  containing  a  solid  and  a 
gas,  the  solid  would  have  to  accept  the  same  amount  of  volumetric  strain  as  the  gas  in  spite  of  the 
fact  that  the  gas  is  much  more  compressible.  The  approximation  correctly  allocates  PdV  work  to 
materials  only  when  compression  occurs  parallel  to  a  material  interface.  This  assumption  is  used 
in  the  MESA,  CTH,  and  KRAKEN  [49]  codes.  Another  approximation  used  for  multiple 
temperatures  and  pressures  in  a  mixed  cell  is  that  all  materials  experience  equal  pressure  change. 
This  approximation  correctly  allocates  PdV  work  among  materials  only  when  a  cell  is  compressed 
perpendicular  to  the  material  interfaces  and  represents  an  opposite  limiting  case  to  constant 
volumetric  strain  approximation.  This  approximation  was  recently  added  to  the  CTH  code  as  a  user 
selectable  option  [50]. 

4.  Conclusions 

Various  factors  related  to  material  interfaces  in  Lagrangian  and  Eulerian  shock  wave 
propagation  codes,  which  can  lead  to  disagreement  between  computations  and  experience,  have  been 
discussed.  For  Lagrangian  codes,  interfaces  are  well  defined;  however,  special  logic  is  needed  to 
prevent  interpenetration  of  two  or  more  bodies.  These  slide-line  algorithms  can  have  user-defined 
parameters  or  options  that,  when  incorrectly  applied,  can  lead  to  considerable  error  in  a  computation. 
For  Eulerian  codes,  interfaces  are  not  well  defined  and  special  logic  is  needed  if  the  interface  is  to 
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be  known  to  greater  than  one  cell  dimension.  Several  means  of  tracking  the  interfaces  and  problems 
associated  with  them  have  been  discussed.  Additionally,  Eulerian  codes  have  difficulty  treating 
sliding  or  material  separation.  Material  properties  in  mixed  cells  are  treated  in  an  ad  hoc  manner. 
Code  users  should  be  aware  of  the  problems  that  can  arise  in  the  treatment  of  material  interfaces 
when  they  occur  so  as  not  to  attribute  them  as  reality  or  new  physics. 
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